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Abstract Sex-determining systems show a 
striking diversity not only among species, but also 
among populations. In reptiles, sex-determination 
is a continuum, from temperature-dependent sex 
determination (TSD) to genetic sex determination (GSD). 
The multi-ocellated racerunner (Eremias multiocellata) 
is reported to be a cryptic ZZ/ZW chromosomal TSD 
species, with male-biased sex ratios at high temperatures 
in two Gansu populations. However, the generality of 
the sex-determining pattern in different populations 
of this species remains unclear. To investigate the 
mode of sex determination in a population of E. 
multiocellata from the desert steppe of Inner Mongolia, 
we first identified sex chromosomes via comparative 
genomic hybridization (CGH). We then conducted a 
thermal manipulation experiment to determine the 
effect of gestation temperature on offspring sex ratios. 
From the CGH studies we found that lizards from 
the Inner Mongolia population possessed ZZ/ZW sex 
chromosomes. However, our thermal manipulation 
experiment showed that gestation temperature did not 
affect the sex ratio of neonates in this population. In 
combination, these results rule out TSD in the Inner 
Mongolia population of E. multiocellata, and suggest that 
there is widespread geographic variation in the sex- 
determining system of this species. 


“Corresponding author: Prof. Xifeng WANG, from Institute of Zoology, 
Chinese Academy of Sciences, Beijing, China, with her research focusing 
on reptilian immunity and temperature sex determination of reptiles. 
E-mail: wangxifeng@ioz.ac.cn 

Received: 1 December 2020 Accepted: 31 May 2021 


Keywords geographic variation, maternal thermal 
environment, reptile, sex chromosome, sex determination, sex 
ratio, sex reversal 


1. Introduction 


Vertebrates exhibit two broad categories of sex-determining 
systems: genotypic sex determination (GSD) where specific 
genetic elements direct the development of the gonads, and 
environmental sex determination (ESD) where external 
cues, such as temperature (temperature-dependent sex 
determination, TSD), determine the sex of an individual. All 
mammals and birds exhibit GSD. A dominant sex gene on 
the Y chromosome region (Sry) initiates male development 
in mammals (Kashimada and Koopman, 2010; Koopman et 
al, 1990). In birds, the dosage effects of double sex and mab- 
3 related transcription factor 1 (Dmrtl) on the Z chromosome 
control sex determination (Ellegren, 2009; Smith et al., 2007; 
Smith et al., 2009). In contrast, sex determination systems are 
highly variable among ectotherms. Both ZZ/ZW and XX/ 
XY sex-chromosome systems occur in all these taxa including 
fish, amphibians and reptiles, as well as TSD in certain reptiles 
and fish (Bachtrog et al, 2014; Capel, 2017; Lambert et al., 2019; 
Nemeshazi et al., 2020; Ospina-Alvarez and Piferrer, 2008; 
Janzen and Krenz, 2004). For example, reptiles have multiple 
sex-determining mechanisms. In all snakes, and some lizards 
and turtles (Ezaz et al, 2006a; Ezaz et al, 2006b; Janzen and 
Krenz, 2004; Schwanz et al., 2013) sex determination involves 
GSD (with ZZ/ZW or XX/XY sex chromosomes), whereas 
in all crocodiles, the tuatara, many turtles and several lizard 
species (Rhen and Schroeder, 2010; Shine et al., 2007) nest 
temperature (in oviparious species) or gestation temperature (in 
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viviparous species) is responsible for determining gonadal sex of 
offspring, 

Environmentally influenced sex determination has been 
found in some ectothermic vertebrates with sex chromosomes, 
in which genes and environment interact to determine the 
phenotype of the sex (Holleley et al, 2015). Temperature can 
override the GSD effect and result in sex reversal in some fishes 
(Baroiller and D’Cotta, 2016; Ospina-Alvarez and Piferrer, 2008), 
amphibians (Alho et al., 2010; Bachtrog et al., 2014; Flament, 2016; 
Lambert et al, 2019; Nemeshazi et al., 2020) and reptiles (Holleley 
et al., 2015; Quinn et al., 2007; Radder et al., 2008; Sarre et al., 
2004; Shine et al, 2002). These species have sex chromosomes but 
exhibit thermal sensitivity in the development of gonadogenesis, 
representing mixed sex-determination systems. For example, 
high incubation temperatures reverse genotypic males (ZZ) 
to phenotypic females in the Australian central bearded 
dragon lizard (Pogona vitticeps), which has GSD with female 
heterogamety (ZW) (Holleley et al, 2015; Quinn et al., 2007). 
Analogously, low temperatures override chromosomal sex 
to generate XX male offspring in a montane scincid lizard 
(Bassiana duperreyi) from south-eastern Australia, which has 
male heterogamety (XY) (Radder et al, 2008; Shine et al, 2002). 
Further work on the interactive effects of sex chromosomes 
and early developmental (embryonic or larval) temperatures on 
offspring sex ratios is needed to elucidate the sex-determining 
systems across a greater diversity of lineages in ectothermic 
vertebrates. 

Interestingly, recent studies have demonstrated that sex- 
determining mechanisms may also differ among populations 
within a species. For example, the association between 
phenotypic sex and the underlying genetic factor (linkage 
group 2) was perfect in a northern Swedish population, but 
weak and variable among families in a southern one in the 
common frog (Rana temporaria) (Rodrigues et al., 2016). Three 
Australian lizards, the snow skink (Niveoscincus ocellatus) (XX/ 
XY system), the bearded dragon (ZZ/ZW system) and another 
skink (Bassiana duperreyi) (XX/XY system), are GSD species 
possessing a mixed sex determination system with temperature- 
induced sex reversal in some populations (Castelli et al., 2020a; 
Cunningham et al, 2017; Dissanayake et al, 2021; Hill et al, 2018; 
Pen et al, 2010). In the snow skink, variation in the magnitude 
of between-year temperature fluctuations might lead to natural 
selection on sex determination across altitudes (Cunningham 
et al, 2017; Pen et al., 2010). Similarly, climatic effects on sex 
determination also occur in the Atlantic silverside, Menidia 
menidia, which exhibits an exceptionally high level of clinal 
variation in sex determination across its geographic range, with 
TSD in southern populations, GSD in northern populations, 
and mixed systems in between (Duffy et al, 2015). However, 
whether such intra-specific divergence in sex-determining 
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systems exists in other lineages remains largely unknown. 

The multi-ocellated racerunner (Eremias multiocellata) is 
a small viviparous lizard (about 65 mm snout-vent length 
SVL) with an extensive distribution ranging from Mongolia, 
Kazakhstan, Kyrgyzstan to Xinjiang, Gansu, and Inner 
Mongolia of China, with the highest known population found 
at an altitude of 2900 m above sea level (asl) (Zhao, 1999). The 
chromosomal system of this species was identified by traditional 
Giemsa staining as 2n = 38, with no obvious sex chromosomes. 
However, a recent study detected female-specific amplified 
signals in two Gansu populations by using comparative 
genomic hybridization (CGH), suggesting this species possesses 
a female-heterogametic sex chromosome system (ZZ/ZW) 
(Wang et al, 2015). Interestingly, gestation temperature has been 
shown to affect offspring sex ratio in two Gansu populations 
at elevations of 1400 m and 2900 m asl, with male-biased sex 
ratios at 35°C and female-biased sex ratios at 25°C (Tang et 
al, 2012; Zhang et al., 2010). In addition, compared to the low- 
elevation population from a relatively warmer climate, the 
high-elevation population of E. multiocellata exhibited a less 
skewed offspring sex ratio (Tang et al., 2012). In the Inner 
Mongolia population of E. multiocellata at a lower elevation of 
1036 m asl, the edge of the species’ altitudinal distribution range, 
extreme environments (ie, low precipitation and extreme high 
temperatures) during reproductive season affect offspring sex 
(Wang et al, 2016). The effect of temperature on offspring sex 
varied between years with different precipitation levels (Wang 
et al, 2016), so a controlled experiment is needed to test if high 
temperature itself (without accompanying aridity) affects 
sex determination in the Inner Mongolia population of E. 
multiocellata. Therefore, here we used CGH to identify the sex 
chromosome system, and conducted temperature manipulation 
experiments to investigate the effects of materal thermal 
(gestation) environment on offspring sex. 


2. Materials and Methods 


2.1. Study site and species The E. multiocellata used in this study 
were collected from Shierliancheng countryside (40°12'17'N, 
111°07'43"E; elevation 1036 m), Jungar Banner, in the desert 
steppe region of Inner Mongolia, China. Shierliancheng is 
characteristic of a cold semi-arid climate with an average 
annual temperature of 6-7°C and annual precipitation of 
300-380 mm (http://www.nmic.gov.cn). This lacertid lizard 
is a viviparous species (adult SVL 54.0-67.3 mm, BM 4.7-8.2 
g) that generally occupies arid or semi-arid regions. Female E. 
multiocellata have a gestation period of two months and produce 
a single litter of 2-5 offspring between July and August (Zhao, 
1999). The selected body temperature and critical thermal 
maximum of E. multiocellata in Inner Mongolia population are 
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(352+0.2)°C and (45.1+0.1)’C, respectively (Li et al. 2017). 


2.2. Identification of sex chromosomes Three male and 
three female adult E. multiocellata were used for sex chromosome 
identification, using the established protocol of CGH (Martinez 
et al., 2008). we prepared both female and male metaphase 
spreads of mitotic chromosomes from peripheral blood 
leukocytes following the published method of Ezaz et al (2005). 
Approximately 100 uL blood was used to set up 2 mL culture 
in RPMI1640 (Gibco) supplemented with 10% fetal bovine 
serum, 100 units/mL penicillin, 0.1 mg/mL streptomycin, 40 
ug/mL phytohaemagglutinin M (PHA M). After incubated at 
30°C in 5% CO, incubators for 72 h, we replaced new medium. 
Eight and six hours before harvesting, 35 mg/mL 5-bromo-2- 
deoxyuridine (BrdU; Sigma) and 75 ng/mL colcemid (Roche) 
were added to each culture. Metaphases were fixed in 3:1 
methanol-acetic acid. Cell suspension was dropped onto wet 
cold slides prechilled at 4°C, and slides were dried on a hotplate 
at 60°C and saved at -80°C. 

Total genomic DNA (gDNA) was extracted from whole 
blood following the phenol-chloroform method, and the 
DNA probe was then labeled with SpectrumGreen-dUTP 
and SpectrumRed-dUTP (Enzo Life Science, Farmingdale, 
NY). Metaphase chromosome slides dried for 2 h at 60°C and 
denatured for 2 min at 60°C in 70% formamide and 2xSSC, 
then dehydrated through a 70%, 80% and 100% ethanol series 
for 1 min each, and air-dried at room temperature until 
hybridization. 

The probe mixture containing 20 mg glycogen, 500 ng 
female probe, and 500 ng male probewere co-precipitated 
with 5 mg of boiled E. macquarii female DNA. Following 
an overnight co-precipitation at -20°C, the probe mix was 
centrifuged at 15 000 g at 4°C for 15 min. The supernatant was 
discarded and the probe DNA pellet was resuspended in 40 
uL of 37°C pre-warmed hybridization buffer and resuspended 
for at least 30 min at 37°C. The hybridization mixture was 
denatured at 70°C for 8 min, immediately placed on ice for 
2 min, and 16 pL of the probe mixture was placed as a single 
drop per slide. For the purposes of this experiment we assumed 
that E. multiocellata have a female ZZ/ZW chromosome system, 
and therefore, non-labeled male gDNA was used to compete 
with DNA probes; in the reciprocal experiment, female gDNA 
was the competitor. After a humid hybridization chamber at 
37°C for 3 days, slides were washed at 55°C in 0.4xSSC, 0.3% 
Tween-20 for 2 min, followed by a second wash in 2xSSC, 0.1% 
Tween-20 for 1 min at room temperature. Slides were immersed 
into a DAPI solution with a concentration of 1 pg/mL for 30 
s. This procedure resulted in the respective hybridization of 
the competing DNA probes, with male and female metaphase 
spreads. Finally, we conducted fluorescene microscopy (Zeiss) to 
detect if there was any sex-specific signal in the spreads. 
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2.3. Thermal treatments A total of 60 gravid female and 18 
adult male E. multiocellata were captured and brought back to the 
laboratory in Beijing in early May 2017. This timing of capture 
allowed us to start treating females prior to the temperature- 
sensitive period (TSP) of offspring sex determination, because 
the TSP is in the middle third of embryonic development in 
reptiles (Wibbels et al., 1991). After recording snout-vent length 
(SVL, +0.01 mm) and body mass (BM, +0.001 g), lizards were 
evenly assigned into three temperature treatments: low (29°C), 
moderate (35°C), and high (38°C). Within each treatment, 
lizards were randomly distributed into three terraria (55.6 
cmx42.5 cmx34.6 cm, lengthxwidthxheight) that were placed 
in temperature-controlled incubators, with each terrarium 
containing 6-7 females and 2 males. Incubators were set at 23°C 
for the low (29°C) and moderate (35°C) temperature treatments, 
and at 29°C for the high temperature treatment (38°C). We then 
used electronic heating pads with temperature control devices 
(YWK-F, Shanghai Medical Instruments, China) attached 
to each terrarium to increase the ambient temperature of 
terraria to 29°C, 35°C, and 38°C for the low, moderate, and high 
temperature treatments respectively, from 08:00 to 18:00 h each 
day. Data loggers (iButton, DS1921, Maxim Integrated Products 
Ltd., USA) were used to monitor the thermal environment 
within each terrarium for the duration of the study. 
Additionally, a 12:12 h (light: dark) photoperiod was provided in 
each terrarium from 07:00 to 19:00 h, using full UV-spectrum 
lamps. Each terrarium contained a sand substrate (5 cm deep) 
and several shelter items. Food (mealworms and crickets dusted 
with extra vitamins and minerals) and water were provided 
ad libitum. We recorded body temperatures of active female 
lizards four times on May 25”, June 1“, June 8” and June 15" 
by inserting a probe of UT325 electronic thermal meter into 
cloacae. 


2.4, Female reproduction and neonate traits Females in the 
experimental terraria were palpated for gravidity every three 
days. Individual females containing soft and oval eggs (ie, late 
developmental-stage embryos) were then transferred into small 
boxes with the same temperature as their original thermal 
treatments to give birth. Then, we checked these boxes every 
day for neonates. Once found, neonates were measured (SVL; 
+0.01 mm) and weighed (BM; +0.001 g). We calculated gestation 
length as the days between female capture and the parturition. 
We obtained a total of 104 neonates, of which 22, 40 and 
42 neonates produced by females from the low, moderate 
and high temperature treatment, respectively. In addition, 10 
out of 22 neonates produced by gravid females from the low 
temperature treatment were stillborn. We identified the sex 
of these 104 neonates by the presence/absence of hemipenes. 
A subset of those neonates (n = 34) did not survive after birth 
and were used to validate our morphological sex identification 
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technique, by inspecting the gonads. Intact gonads were 
dissected out, snap-frozen in liquid nitrogen and embedded in 
an optimal cutting temperature compound (OCT) at -20°C 
to -30°C. Frozen gonads were cut into 6 ym slices with the 
microtome portion of the cryostat. Each slice was then placed 
on a glass slide and stained with hematoxylin and eosin (H&E 
stain). Sex was identified under an optical microscope based on 
standard testis and ovarian morphology (Figure 1). The sex of 
all of neonates identified by gonadal inspection was consistent 
with our morphological sex-identification technique (everting 
hemipenes; Figure 1); subsequently, neonate sex was identified 
using presence/absence of hemipenes. 


2.5. Statistical analysis We used linear models (LM; “stats” 
package) (R Core Team., 2019) to detect the differences of 
maternal SVL and BM, and the ambient temperatures 


Female 
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inside terraria between thermal treatments. We used the 
generalized linear models (GLM; “stats” package) (R Core 
Team., 2019) to detect the effects of thermal treatments on 
female reproductive success using a binomial distribution with 
a logit link. Generalized least squares models (GLS; “nlme” 
package) (Pinheiro et al, 2017) were used to detect the effects 
of thermal treatments on gestation period and litter mass 
with maternal SVL as a covariate, because within-group 
variance differed between treatment groups, we allowed 
for treatment-dependent variances in the GLS model by the 
varldent function. The Kruskal-Wallis test was used to detect 
the effects of thermal treatments on litter size (the total number 
of offspring born by a female including stillborn). We used 
linear mixed-effects models (LME “nlme” package) (Pinheiro 
et al., 2017) to test the effects of thermal treatments on neonate 
SVL and BM with maternal ID as a random factor. The same 


Male 


Figure 1 Methods of sex identification for Eremias multicellata used in this study. (A) Microscopy of frozen sections of gonads from E. 
multicellata. (B) Manual eversion of hemipenes. Sex of E. multicellata was determined by presence/absence of hemipenes. Arrow points to 


the hemipenes. 
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model LME was used to detect the effects of thermal treatments 
on body temperature of the E. multiocellata. In order to detect 
the effects of thermal treatments on neonate survival, total 
neonate (including both surviving and dead individuals) sex 
ratio, sex ratio of neonates that were dead after birth, and sex 
ratio of neonates that survived after birth, we used generalized 
linear mixed effect models (GLMMs; “Ime4” package) (Bates 
et al, 2015) with a binomial distribution and a logit link, and 
maternal ID as a random factor. We compared each model 
including thermal treatments with a model that only contained 
the random factor using the “anova” function to get an overall 
chi-squared and P-value for the thermal treatments effect. Post- 
hoc tests for the analysis of litter mass and gestation period 
were performed by emmeans function with “tukey” method 
(“emmeans” package) (Lenth R.V., 2018) and we used pairwise 
Mann-Whitney tests to carry out post-hoc tests for significant 
Kruskal-Wallis test and corrected the P-values by the “fdr” 
method (“stats” package) (R Core Team., 2019). We checked 
the GLMMs with binomial distribution for the dispersion 
parameter by “dispersion_glmer’ function (“blmeco’ package) 
(Korner-Nievergelt et al., 2015), and the dispersion parameters of 
all binomial models are (close to) 1. We conducted all statistical 
analyses in R version 3.5.1 (R Core Team, 2019). 


3. Results 


3.1. Sex chromosomes of E. multiocellata We detected a 
significantly enhanced and specific hybridization signal in 


SpG-male probe 
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one of the microchromosomes in cells from all females, but 
not in males (Figure 2). This female-specific chromosome is 
therefore, by definition, a W chromosome, which means that E. 
multiocellata from the Inner Mongolia population possess a ZZ/ 
ZW sex chromosome system. 


3.2. Female body temperature and reproduction The 
ambient temperatures inside terraria were significantly 
different among the three temperature treatments (F, s1 = 
126.59, P < 0.001; Figure 3A). Correspondingly, daily average 
body temperatures of gravid females from 08:00 to 18:00 
were also significantly different among thermal treatments 
(F, 36 = 247.4, P < 0.001). Mean (+SE) body temperatures of 
gravid females were as follows: high temperature treatment 
(37.9+40.1°)C, moderate temperature treatment (34.6+0.1)’C, low 
temperature treatment (28.9+0.1)°C (Figure 3B). 

Gestation period of gravid females was significantly 
affected by thermal treatment, decreasing as female gestation 
temperature increased (Table 1 and Figure 4). Female 
reproductive success (the proportion of females producing live 
offspring) was 65.0% (13/20) for females in the high temperature 
group, 80.0% (16/20) in the moderate temperature group, and 
50.0% (10/20) in the low temperature treatment (Figure 4). 
This difference in reproductive success among temperature 
treatments was not statistically significant (Table 1). Thermal 
treatment significantly affected litter mass and litter size. 
Females from the high temperature treatment had greater litter 
mass and larger litter size than those from the low temperature 
treatment (Table 1 and Figure 4). 


SpG-female probe Merged 


Female spread 


Male spread 


Figure 2 Results of comparative genomic hybridization (CGH) in female (upper) and male (lower) metaphases from Eremias multiocellata 
in our study. (A, E) DAPI stained metaphase chromosome spread; (B, F) Metaphase hybridized with SpectrumGreen-dUTP labeled male 
whole genomic DNA; (C, G) Metaphase hybridized with SpectrumRed-dUTP labeled female whole genomic DNA; (D, H) Merged images 
of metaphase hybridized with male and female whole genomic DNA. Scale bar indicates 100 pm. Arrow indicates W chromosome. 
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3.3. Neonate phenotypes and sex ratio Neither neonate SVL 
nor BM differed among thermal treatments (Table 1 and Figure 
5). The proportion of neonates that did not survive after birth 
was higher in the high (78.6%) and low (75.0%) temperature 
treatment than the moderate (45.0%) temperature treatments 
(Figure 5). 

Most notably, female gestation temperature did not affect 
neonate sex ratio (Table 1). The male ratio of neonates (54.5%, 
52.5% and 52.4% for the low, moderate and high temperature 
treatments, respectively, which was very close to 50:50 in all 
three groups, did not differ between the three temperature 
treatments (Figure 5 and Table 1). 


4. Discussion 


Our study revealed that E. multiocellata from the Inner Mongolia 
population have a female ZZ/ZW chromosome system, and 
that female gestation temperature does not affect offspring 
sex ratio; this is inconsistent with the TSD pattern found in 
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Figure 3 Experimental thermal treatments. (A) Ambient 
temperatures recorded in the terraria from the three 
temperature treatments, and (B) active body temperatures 
recorded during gestation in gravid female Eremias multicellata. 
Values are expressed as mean+SE. 
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the two Gansu populations of E. multiocellata (Tang et al, 2012; 
Zhang et al., 2010). This between-population discrepancy in 
sex-determining system adds another example of intraspecific 
variation in TSD and GSD patterns in reptiles, in addition to the 
phenomenon first reported in a viviparous Australian montane 
lizard, the snow skink, Niveoscincus ocellatus (Cunningham et 
al, 2017; Pen et al., 2010). Sex-determining systems are quite 
variable in reptiles, involving environmental and genetic sex 
determination, and only some reptiles have been identified to 
possess sex chromosome differentiation (Bachtrog et al., 2014; 
Georges et al., 2010; Martinez-Juarez and Moreno-Mendoza, 
2019). Traditional Giemsa staining methods did not elucidate 
sex chromosomes in E. multiocellata because the traditional way 
identifies chromosome pairs based on morphological differences 
(Wang et al. 2015). Instead, we used CGH, which distinguishes 
sex chromosomes through the non-homologous region, to 
reveal that E. multiocellata from the Inner Mongolia population 
had ZZ/ZW sex chromosomes (Figure 2). The sex chromosomes 
of some reptiles are extremely cryptic, and can only be detected 
via such high-resolution cytogenetic techniques (Badenhorst 
et al, 2013; Ezaz et al, 2005; Ezaz et al., 2006b; Kawai et al., 2007; 
Martinez et al, 2008). The Gansu populations of E. multiocellata 
also have ZZ/ZW sex chromosomes, but the sex was affected 
by temperature (Wang et al., 2015). These results suggest that 
the E. multiocellata populations we have studied all have sex 
chromosomes, but that the genetic effects on offspring sex 
can be overridden by extreme temperatures in the Gansu 
populations (Tang et al, 2012; Zhang et al, 2010). 

The gestation temperatures we used in our experiment, 
which mimicked female body temperatures measured in the 
field during the reproductive season (Li et al, 2017), did not 
affect offspring sex ratio in the Inner Mongolia population 
of E. multiocellata. However, a field warming experiment 
on this population showed that females produced male- 
skewed offspring under extremely high temperatures and 
dry environment, but not in response to extremely high 
temperatures alone (Wang et al., 2016). The sex-ratio bias 
found under warm and dry conditions could have been due 
to water deprivation or maternal manipulation of offspring 
sex ratios in response to nutrient deficiency caused by the 
extreme conditions. For example, water deprivation during 
pregnancy could induce a male-biased secondary sex ratio 
in Aspic viper and Zootoca vivipara (Dupoué et al., 2019). In 
addition, in an agamid lizard (Amphibolurus muricatus) from 
Australia, females fed a poor-quality diet produced highly male- 
biased sex ratios (Warner and Shine, 2007). These instances 
indicate that female-to-male sex reversal may be triggered 
by stressors, perhaps due to the effects of glucocorticoid 
hormones and reactive oxygen species on sex determination 
(Castelli et al, 2020b). In addition, the ecological driving forces 
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Table 1 Differences between gestation temperature treatments in maternal and neonate traits in the multi-ocellated racerunner Eremias 
multicellata from Inner Mongolia. 


. Treatment wipe : 
Variable 29°C 35°C 38°C Statistical analysis 
Maternal SVL (mm) 65.108+1.017 64.468+0.774 65.068+1.147 F135 = 0.147, P = 0.864 
Maternal body mass (g) 5.615+0.179 5.576+0.201 5.890+0.332 F 335 = 0.472, P = 0.628 
F, 3 = 119.808, P < 0.001 
Gestation period (day) 60.100+2.681 45.560+1.494 28.230+0.533 29°C*>35°C°>38°C° 
a-b: P = 0.001; a-c: P < 0.001; b-c: P< 0.001 
Reproductive success 50.0% (10/20) 80.0% (16/20) 65.0% (13/20) X = 4.054, df = 2, P = 0.132 
Fy = 17.885, P < 0.001 
Litter mass (g) 1.021+0.031 1.292+0.070 1.755+0.159 29°C*<35°C'<38°C° 
a-b: P = 0.004; a-c: P = 0.001; b-c: P = 0.032 
xX = 9.131, df=2, P = 0.010 
Litter size 2.200+0.133 2.500+0.158 3.230+0.281 29°C* <38°C 
a-b: P = 0.217; a-c: P< 0.019; b-c: P = 0.056 
Neonate SVL (mm) 29.319+0.247 29.257+0.208 29.172+0.204 F, 36 = 0.201, P= 0.819 
Neonate body mass (g) 0.498+0.011 0.51740.013 0.543+0.012 F, 36 = 1.589, P = 0.218 
Neonate survival ratio 75.0% (9/12) 45.0% (18/40) 78.6% (33/42) X = 2.240, df = 2, P = 0.326 
Neonate sex ratio (3%) 54.5%(12/22) 52.5% (21/40) 52.4% (22/42) xX = 0.029, df = 2, P =0.986 
Neonate sex ratio 6 2 
(Deadarecebath oo) 66.7% (2/3) 63.6% (14/22) 55.6% (5/9) IF a U267, aif = 2 PS OLID 
Neonate sex fatio 55.6% (5/9) 38.9% (7/18) 51.5% (17/33) X =0.971, df = 2, P = 0.615 


(Survived after birth, 4%) 


The maternal body size (SVL and BM) were analyzed with linear models; the female reproductive success were analyzed with generalized 
linear models; the gestation period and the litter mass were analyzed with linear generalized least squares models; the litter size were 
analyzed with Kruskal-Wallis test; the neonate SVL and BM were analyzed by LME with mother ID as random factor; neonate survival and 
neonate sex ratio were analyzed with generalized linear mixed effect models. 
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Figure 4 Differences between gestation temperature treatments in maternal snout-vent length (SVL), body mass (BM), gestation period 
length, reproductive success, litter size and litter mass in the multi-ocellated racerunner Eremias multicellata from Inner Mongolia. 
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Figure 5 Neonate snout-vent length (SVL), body mass (BM), and sex ratio in different thermal treatments during maternal gestation of 
the multi-ocellated racerunner Eremias multicellata from Inner Mongolia. 


of different sex determining systems among populations of E. 
multiocellata remain unrevealed, but could be related to among- 
population differences in climate conditions. Our study site at 
Inner Mongolian has lower ambient temperatures in terms 
of average, the lowest and highest temperatures and higher 
precipitation than two sites of Mingin and Tianzhu in Gansu 
Province (Table 1). The dry and warm condition in Mingin and 
Tianzhu might be a reason for the evolution of a mixed sex- 
determining system with temperature-induced sex reversal in 
these populations. 

Currently there is a knowledge gap on intraspecific 
variation in sex-determining systems, concerning the 
evolutionary driver of different sex-determination patterns 
among populations. For example, climate or elevation did 
not explain within-species variation in sex reversal in the 
Australian dragon lizard (Cornejo-Paramo et al., 2020a) or sex- 
chromosome differentiation in common frogs (Phillips et al., 
2020). The population divergence of sex-determining systems 
found in N. ocellatus was attributable to the increased rate of 
female maturation in the lowland population, and the higher 
magnitude of between-year temperature fluctuations in the 
highland population (Pen et al., 2010). Therefore, climatic 


effects on lizard life history may induce divergent natural 
selection forces on sex determination across altitudinal clines 
eg, a relatively long activity season favours an evolutionary 
shift from GSD to TSD in lowland populations of N. ocellatus 
(Cunningham et al., 2017; Pen et al., 2010). However, this 
explanation cannot be applied to the E. multiocellata system, 
because TSD occurs in high-altitude populations rather than 
low-altitude populations in this species. In addition, a meta- 
analysis on sex-determining systems of non-avian reptiles 
falsified the assumption that populations living in either 
highly variable or cold climatic conditions should evolve GSD 
to buffer the populations from extreme sex ratios (Cornejo- 
Páramo et al., 2020b). The selective forces acting on the evolution 
of TSD-GSD systems remain a mystery and warrants further 


study in future. 
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Appendix 


File 1 Excel 1 presented the data of female gestation traits including maternal SVL, BM, Gestation period, litter size, litter mass and 
reproductive success in the multi-ocellated racerunner Eremias multicellata from Inner Mongolia with different temperature treatment. Excel 
2 presented the data of neonate traits (SVL, BM, sex and survival) with different temperature treatment of their mother. The file can be 
downloaded from the website https://pan.baidu.com/s/1M8MuOyWLkyzBH_rp_S9sqw (access code: k9og). 


File 2 Meteorological data of Tianzhu, Mingin and Eerduosi (Inner Mongolia) from http://data.cma.cn in recent 20 years. LowestT means 
the lowest temperature, highestT represents the highest temperature, and averageT means the average temperature. RH represents relative 
humidity. The file can be downloaded from the website https://pan.baidu.com/s/1AZ30W91K-lij3LHjXsZ1BQ (access code: v9n5). 


